High-energy 7-rays from stellar associations 
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ABSTRACT 

It is proposed that TeV 7-rays and neutrinos can be produced by cosmic rays 
(CRs) through hadronic interactions in the innermost parts of the winds of mas- 



CD 
(N 

O 

Q 



> 
00 

(N 

m 
o 



sive and B stars. Convection prevents low-energy particles from penetrating 
into the wind, leading to an absence of MeV-GeV counterparts. It is argued that 
groups of stars located close to the CR acceleration sites in OB stellar associations 
may be detectable by ground-based Cerenkov telescopes. 

Subject headings: gamma rays: observations — gamma rays: theory — stars: early- 
type 



Introduction 



Several 7-ray sources are thought to be related with early-type stars and their neigh- 
borhoods (e.g., Montmerle 1979; Casse & Paul 1980; Bykov & Fleishman 1992a,b; Bykov 
2001, Romero & Torres 2003). Recently, the first (and only) TeV unidentified source was 
detected in the Cygnus region (Aharonian et al. 2002), where a nearby EGRET source (3EG 
J2033+4118) has a likely stellar origin (White & Chen 1992; Chen et al. 1996; Romero et 
al. 1999; Benaglia et al. 2001). Here, we explore whether CR illumination of stellar winds 
of O and B stars can lead to Galactic TeV 7-ray sources. 

2. The model 

O and B stars lose a significant fraction of their mass in stellar winds with terminal 
velocities V^o ~ 10 3 km s _1 . With mass loss rates as high as M* = (10 -6 — 10~ 4 ) M & 
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yr -1 , the density at the base of the wind can reach 10~ 12 g cm -3 (e.g., Lamers & Cassinelli 
1999, Ch. 2). Such winds are permeated by significant magnetic fields, and provide a matter 
field dense enough as to produce hadronic 7-ray emission when pervaded by relativistic 
particles. A typical wind configuration (Castor, McCray, & Weaver 1975; Volk & Forman 
1982; Lamers & Cassinelli 1999, Ch. 12) contains an inner region in free expansion (zone I) 
and a much larger hot compressed wind (zone II). These are finally surrounded by a thin 
layer of dense swept-up gas (zone III); the final interface with the interstellar medium (ISM). 
The innermost region size is fixed by requiring that at the end of the free expansion phase 
(about 100 years after the wind turns on) the swept-up material is comparable to the mass 
in the driven wave from the wind, which happens at a radius -R w ind = Ko(3M^/47rp Kx)) 1 ^ 2 ) 
where po ~ m p no is the ISM mass density, with m p the mass of the proton and no the 
ISM number density. After hundreds of thousands of years, the wind produces a bubble 
with a radius of the order of tens of parsecs, with a density lower (except that in zone I) 
than in the ISM. In what follows, we consider the hadronic production of 7-rays in zone I, 
the innermost and densest region of the wind. The matter there will be described through 
the continuity equation: M+ = 4nr 2 p(r)V (r) , where p(r) is the density of the wind and 
V(r) = Voo(l — Ro/r) 13 is its velocity. is the terminal wind velocity, and the parameter (3 
is ~ 1 for massive stars (Lamers & Cassinelli 1999, Ch. 2). R is given in terms of the wind 
velocity close to the star, V ~ lO -2 !/^, as R = R+(l — (K)/V 00 ) 1 / /3 ). Hence, the particle 
density is n(r) = M*(l — R / 'r)" 13 / ' (AixmpV^r 2 ) . 

Not all CRs will enter into the base of the wind. Although wind modulation has 
only been studied in detail for the case of the relatively weak solar wind (e.g. Parker 
1958, Parker & Jokipii 1970, Kota & Jokipii 1983, Jokipii et al. 1993), a first approach 
to determine whether particles can pervade the wind is to compute the ratio (e) between 
the diffusion and convection timescales: = 3r 2 /D and t c = 3r/V(r), where D is the 
diffusion coefficient, and r is the position in the wind. Only particles for which e < 1 
will be able to overcome convection and enter the dense wind region to produce 7-rays 
through pp interactions. The diffusion coefficient is D ~ A r c/3, where A r is the mean- 
free-path for diffusion in the radial direction. As in White (1985) and Volk and Forman 
(1982), the mean-free-path for scattering parallel to the magnetic field (B) direction is 
assumed as Ay ~ I0r g = lOE/eB, where r g is the particle gyro-radius and E its en- 
ergy. In the perpendicular direction A is shorter, Aj_ ~ r g . The mean-free-path in the 
radial direction is then given by A r = Aj_ 2 sin 2 9 + Ay 2 cos 2 9 = r g (10 cos 2 9 + sin 2 9), where 
cos -2 9 = 1 + (B^/Br) 2 . Here, the geometry of the magnetic field is represented by the 
magnetic rotator theory (Weber and Davis 1967; White 1985; Lamers and Cassinelli 1999, 
Ch. 9) B^/Br = (T4/V r 0O )(l + r/i^) and B r = B*(Rjr) 2 , where V* is the rotational velocity 
at the surface of the star, and B* the surface magnetic field. Using all previous formulae, e ~ 
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Fig. 1. — Left: Minimum proton energy needed to overcome the wind convection at different 
distances from the star. Here V* = 250 km s -1 , = 1750 km s _1 , -R* = 12 R Q . Right: 
Opacities to pair production as a function of the 7-ray energy for different creation places 
R c . Here, the star has T e s = 38000 K. 

{ZeB+V^r - R*){Rjr)\l + (VJV^I + r/i?,)) 2 ) 3 / 2 } / {E p c(10 + (VjV^l + r/i? + )) 2 )} . The 
latter equation defines a minimum energy E™ m (r) below which the particles are convected 
away from the wind (shown in Fig. 1, left panel). Note that E™ m (r) is an increasing function 
of r, so particles that are not convected away in the outer regions of the wind are able to 
diffuse up to its base. £™ in (r) can then be effectively approximated by E™ ul (r ^> /?*) in 
subsequent computations. Only particles with energies higher than a few TeV will interact 
with nuclei in the inner wind and ultimately generate 7-rays, substantially reducing the flux 
in the MeV - GeV band. 

The opacity to pair production of the 7-rays in the UV stellar photon field can be 
computed as t(R c ^E 1 ) = J °° N(E ic )a e - e +(E ir ,E 1 )dE iK dr ) where E± is the energy of the 
photons emitted by the star, E 1 is the energy of the 7-ray, R c is the place where the photon 
was created within the wind, and a e - e +(E±, i? 7 ) is the cross section for 77 pair production 
(Cox 1999, p. 214). The stellar photon distribution is that of a blackbody peaking at typical 
star effective temperatures (T cS ), N(E+) = (irB(E*)/hE*c)Rl/r 2 , where h is the Planck 
constant, and B{E±) = (2E* / (he) 2 ) / {e E *l kT - s - 1). t(R c ,E^) is shown in Fig. 1 (right 
panel) for different photon creation sites (R c <^ R W md)- 7-ray photons of TeV and higher 
energies do not encounter significant opacities in their way out of the wind, unless they 
are created at its very base, hovering over the star (which is unlikely to happen because 
-Rwmd R± and the proton propagates in a high magnetic field environment). Although 
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we show the opacity for values of the photon energy as low as 100 GeV, most of the 7- 
rays will have higher energies, since only protons with E p > E™ m will enter the wind. The 
grey (light-grey) box in the figure shows typical energies of 7-rays for the case of a surface 
magnetic field B± = 10 G (100 G). There is a large uncertainty about the typical values for 
the magnetic field in the surface of O and B stars, but recent measurements favor B± > 100 
G (e.g., Donati et al. 2001; 2002). 



3. 7-ray and neutrino emission 



The differential 7-ray emissivity from 7r°-decays can be approximated by q 1 (E 1 ) = 
A7ia pp (E p )(2Z ( p a l 7T0 /a) J p (E y )i]xQ(E p — E' p mn ) at the energies of interest. The parameter 
i]a takes into account the contribution from different nuclei in the wind (for a standard 
composition t/a ~ 1-5, Dermer 1986). J P (E y ) is the proton flux distribution evaluated at 
E = E 7 (units of protons per unit time, solid angle, energy-band, and area). The cross 
section a pp (E p ) for pp interactions at energy E p « 10E^ can be represented above E p ps 10 
GeV by app(Ep) « 30 x [0.95 + 0.06 log(£ p /GeV)] mb (e.g., Aharonian & Atoyan 1996). 
Z^l^o is the so-called spectrum-weighted moment of the inclusive cross-section. Its value for 
different spectral indices a is given, for instance, by Drury et al. (1994). Finally Q(E P — E^ m ) 
is a Heaviside function that takes into account the fact that only CRs with energies higher 
than i?™ m (r ^> R+) will diffuse into the wind. The spectral 7-ray intensity (photons per unit 
time and energy-band) is Ij(E 7 ) = / n(r)q^(E^)dV, where V is the interaction volume. The 
luminosity in a given band is L 7 = J^ ind /Jj 2 n(r) q J (E^)E 1 (4irr 2 )dr dE y (e.g. Torres et al. 
2003, Romero et al. 2003 for details). Assuming a canonical spectrum for the relativistic CR 
population, J P {E P ) = (c/4tc)N(E p ) = (c/An)K p E p ~ a , the result (in the range E 1 ~ 1 — 20 
TeV) can be expressed in terms of the normalization K p and will depend on all other model 
parameters, mainly on the proton (photon) spectral index, the ISM density, the terminal 
velocity, and the mass-loss rate. Very mild dependencies appear with (3 and Table 
1 presents results for the luminosity for typical values of all these parameters. We have 
fixed = 1O~ 5 M yr _1 , /3 — 1, and i? + = 12i? Q in this example. The mass contained 
in the innermost region of the wind, M win d, is also shown. L 7 ~ 10 25 ~ 30 K P erg s _1 can be 
obtained as the luminosity produced by one particular star; the total luminosity of a group of 
stars should add contributions from all illuminated winds. Convolving the previous integral 
with the probability of escape (obtained through the opacity as e~ r ) does not noticeably 
change these results. Finally, it is possible to factor out the normalization in favor of the 
CR enhancement in the region where the wind is immersed. The CR energy density is 
^cr — J N p (E p )E p dE p = 9.9K P eV cm -3 = <^cr,o, where q is the enhancement factor of 
the CR energy density with respect to the local value, wcr,o (energies between 1 GeV and 
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Table 1: Examples for hadronic 7-ray luminosities from typical stellar wind configurations. 



Model 




n 




-^^wind 


L« =1 - 9 /K p 


L a = 2 -°/K p 






(km s- 1 ) 


(cm -3 ) 


(pc) 


(M e ) 


(erg s _1 ) 


(erg s _1 ) 


(erg s _1 ) 


a 


1750 


10 


0.07 


0.0004 


2 xlO 28 


7 xlO 26 


3 xlO 25 


b 




1 


0.24 


0.0013 


5 xlO 28 


2 xlO 27 


8 xlO 25 


c 




0.1 


0.75 


0.0041 


2 xlO 29 


7 xlO 27 


3 xlO 26 


d 




0.01 


2.4 


0.0130 


5 xlO 29 


2 xlO 28 


8 xlO 26 


e 


1000 


10 


0.09 


0.0009 


4 xlO 28 


1 xlO 27 


6 xlO 25 


f 




1 


0.31 


0.0030 


1 xlO 29 


5 xlO 27 


2 xlO 26 


g 




0.1 


0.99 


0.0095 


4 xlO 29 


1 xlO 28 


6 xlO 26 


h 




0.01 


3.1 


0.0301 


1 xlO 30 


5 xlO 28 


2 xlO 27 


i 


800 


10 


0.11 


0.0013 


5 xlO 28 


2 xlO 27 


9 xlO 25 


J 




1 


0.35 


0.0042 


1 xlO 29 


7 xlO 27 


3 xlO 26 


k 




0.1 


1.1 


0.0133 


5 xlO 29 


2 xlO 28 


9 xlO 26 


1 




0.01 


3.5 


0.0421 


1 xlO 30 


7 xlO 28 


3 xlO 27 



20 TeV). Then, K p ~ (0.2 - 0.3)?. 

The i^ + P^ neutrino flux (F V (E V )) will be derived from the observed 7-ray flux (F 7 (£" 7 )) 
by imposing energy conservation (see Alvarez-Muniz & Halzen 2002 for details): J E^F^E^dE^ = 
C f E u F u (E u )dE u , where the limits of the integrals are E™ 1 ^ (E™^), the minimum (max- 
imum) energy of the photons [neutrinos] and the pre-factor C is a numerical constant of 
order one. Using the resulting v flux, the signal for the detection of //-events can be 
approximated as (Anchordoqui et al. 2003) S = T obs J dE l/ A eS F u (E I/ )P u ^ fl (E u ) whereas 
the noise will be given by iV = (T obs / dE u A cS F B {E u )P u ^(E u )AQ) (1/2) , whereAfi is the 
solid angle of the search bin (A^i°xi° « 3 x 10' 4 sr for ICECUBE, Karle 2002) and 
Fb(E v ) < 0.2 (E u /GeV)~ 3 - 21 GeV -1 cm' 2 s _1 sr -1 is the i^ + P^ atmospheric z/-flux (Volkova 
1980, Lipari 1993). Here, P V ^^(E V ) « 3.3 x 10" 13 (E u /GeV) 2 - 2 denotes the probability that a 
v of energy E v ~ 1 — 10 3 GeV, on a trajectory through the detector, produces a muon (Gaisser 
et al. 1995). T Q b s is the observing time and A e g the effective area of the detector. Those 
systems producing a detectable 7-ray flux above 1 TeV are prime candidates to also be 
detectable neutrino sources, see below. 
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4. Source location and luminosity 

The flux expected at Earth from an isolated star can be computed as F 1 (E 1 > ITeV) = 
(1/47T-D 2 ) J^ wind J 1TeV n(r) qj(Ej) Attt 2 dr dEj. The models in Table 1, at 2 kpc, give fluxes 
in the range (1 x 10~ 20 — 7 x 10~ 16 )i^ p photons cm -2 s _1 . Hence, there are models for which 
a small group of ~ 10 stars in a region with a CR enhancement factor of ~ 100 might be 
detectable at the level of ~ 10~ 14 photons cm -2 s -1 . 

CRs are expected to be accelerated in OB associations through turbulent motions and 
collective effects of stellar winds (e.g. Bykov & Fleishman 1992,b). The main acceleration 
region for TeV particles would be in the outer boundary of the supperbubble produced 
by the core of the association. If there is a subgroup of stars located at the acceleration 
region, their winds might be illuminated by the locally accelerated protons, which would 
have a distribution with a slope close to the canonical value, a ~ 2. For stars out of 
the acceleration region, the changes introduced in the proton distribution by the diffusion 
of the particles (a steepening of its spectrum) would render the mechanism for TeV 7-ray 
production inefficient. This can be seen from Table 1 through the strong dependency of the 
predicted TeV luminosity on the spectral slope of the particles. 

An important assumption in our model is that the diffusion coefficient is a linear function 
of the particle energy in the inner wind. This is indeed an assumption also in both Volk 
& Forman (1982) and White's (1985) models of the particle diffusion in the strong winds 
of early-type stars, among other studies. Measurements of the solar wind, however, seem 
to suggest a harder relation with energy (e.g., D oc E 0A ~ - 5 , Ginzburg & Syrovatskii 1964, 
p. 336). If such a relation would hold for the inner wind of an O star (where pp interactions 
occur), depending on the constant of proportionality, it could yield a higher value of E™ in 
and hence a lower 7-ray luminosity. However, contrary to what happens with the Sun, in 
early-type stars line-driven instabilities are expected to produce strong shocks in the inner 
wind (Lamers & Cassinelli 1999). In such a scenario, as emphasized by White (1985), to 
expect that particles will diffuse according to the Bohm parameterization seems not to be 
unreasonable. As we discuss in the next section, direct observation of TeV sources of stellar 
origin can shed light on the issue. 

5. Application: the unidentified TeV source 

The HEGRA detection in the vicinity of Cygnus OB2, TeV J2032+4131 (Aharonian 
et al. 2002), presents an integral flux F 7 (£ 7 > ITeV) = 4.5(±1.3) x 10~ 13 photons 
cm -2 s _1 , and a 7-ray spectrum F 7 (£' 7 ) = B(E 1 /TeV)~ r photons cm -2 s _1 TeV -1 , where 
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B = 4.7(±2.1 stat ± 1.3 sys ) x 10~ 13 and T = 1.9(±0.3 sta t ± 0.3 sys ). No counterparts at lower 
energies are presently identified (Butt et al. 2003, Mukherjee et al. 2003). The source flux 
was constant during the three years of data collection. The extension of the source (5.6±1.7 
arcmin) disfavors a pulsar or active galactic nuclei origin. The absence of an X-ray coun- 
terpart additionally disfavors a microquasar origin. Instead, the location of the TeV source, 
separate from the core of the association, and coincident with a significant enhancement of 
the star number density (see Fig. 1 of Butt et al. 2003) might suggest the scenario outlined 
in the previous section. 

A nearby EGRET source is, on the other hand, coincident with the center of the asso- 
ciation, where it might be produced either in the terminal shocks of powerful stars therein 
existing (White and Chen 1992, Chen et al. 1996), or in the colliding wind binary system 
Cyg OB2 #5 (Benaglia et al. 2001), or in a combination of these scenarios. Contributions 
from the inner winds of OB stars as in the model herein explored cannot be ruled out. These, 
however, are not expected to dominate because of wind modulation (at low energies) and of 
the softening of the CR spectrum while diffusing from the superbubble accelerating region, 
which significantly diminish the number of pp interactions in the winds. 

Our model could explain the unidentified TeV source without requirements other than 
the presence of the already observed stars and a reasonable CR enhancement if the density 
of the original ISM was rather low. Butt et al. (2003) argued for a density of no = 30 
cm -3 . However, this should be taken as a generous upper limit: a) Apparently, there is 
no star formation currently active at the position of the source, b) The CO-H 2 conversion 
factor used to compute the density has been taken as the normal Galactic one, but it could 
be lower in the neighborhood of star forming environments (e.g., Yao et al. 2003). c) 
The particle density within the TeV source region has been averaged from a velocity range 
integrated along the line of sight corresponding to 3700 pc and including the core of the 
Cygnus association, d) The TeV source will actually be immersed in the zone II of the winds 
of the several powerful stars therein detected, which should have swept the ISM away and 
diminished its density. Our models (e.g., model g of Table 1), which in fact take for the 
stellar parameters an average value from the stars in Table 3 of Butt et al. (2003), show that 
the illumination of the innermost regions of the winds of ~ 10 stars with a CR enhancement 
of ~ 300 in a medium density of about 0.1 cm -3 may be enough to produce the HEGRA 
detection. The neutrino flux that results from a hadronic production of the TeV 7-ray source 
would not produce a significant detection in AMANDA II, which is consistent with the latest 
reports by the AMANDA collaboration (Ahrens et al. 2003). In ICECUBE, however, the 
signal-to-noise is ~ 1.8 for 1 yr of observation (for energies above 1 TeV, an effective area 
of 1 km 2 , before taking into account neutrino oscillations effects). If ICECUBE can reach a 
1° x 1° or finer search bin, and a km 2 effective area at TeV energies, a long integrating time 
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could distinguish the hadronic origin of the HEGRA detection. 

6. Concluding Remarks 

Hadronic interactions within the innermost region of the winds of O and B stars can 
produce significant 7-ray luminosities at TeV energies, with low brightness at other energies. 
At distances less than a few kpc, several illuminated winds pertaining to subgroups of stars 
located at CR acceleration regions in OB associations might be detected by Cerenkov tele- 
scopes. A reasonable set of model parameters can be found to produce a flux compatible 
with the only unidentified TeV source known. A candidate selection for possible new TeV 
sources, based on these predictions, will be reported elsewhere. 
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